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Abstract 

We compare dewetting characteristics of a thin nonwetting sohd film in the absence of stress, 
for two models of a wetting potential: the exponential and the algebraic. The exponential model 
is a one-parameter (r) model, and the algebraic model is a two-parameter (r, m) model, where 
r is the ratio of the characteristic wetting length to the height of the unperturbed film, and m is 
the exponent of h (film height) in a smooth function that interpolates the system's surface energy 
above and below the film-substrate interface at z = 0. The exponential model gives monotonically 
decreasing (with h) wetting chemical potential, while this dependence is monotonic only for the 
m = 1 case of the algebraic model. Linear stability analysis of the planar equilibrium surface is 
performed. Simulations of the surface dynamics in the strongly nonlinear regime (large deviations 
from the planar equilibrium) and for large surface energy anisotropies demonstrate that for any m 
the film is less prone to dewetting when it is governed by the algebraic model. Quasiequilibrium 
states similar to the one found in the exponential model exist in the algebraic model as well, 
and the film morphologies are similar. 

PACS numbers: 68.55.-a 
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I. INTRODUCTION 



Dewetting of lattice-matched ultrathin solid films (such as the sub-10 nm Si film on the 
Si02 substrate) was recently observed in experiments at temperatures around 800°C 
Presumably, the cause for film dewetting is a long-range, attractive film-substrate interaction 
(also called wetting interaction) which amplifies perturbations of the planar film surface and 
makes the film height decrease locally until the surface reaches the substrate, resulting in 
the formation of an array of islands. At this most general level of description dewetting of 
solid films is similar to dewetting of liquid films (which has been studied for many years 
3|, |j|), the only difference is the nature of the mass transport, i.e. the thermally activated 
surface diffusion of adatoms in the former case vs. the fluid flow in the latter case. There is, 
however, two determinative reasons of as to why the dynamics of dewetting in these systems 
is qualitatively different. One reason is the nonzero (and generally, strong) anisotropy of the 
solid film surface energy (tension) which is not present in liquids. As has been shown by the 
author in Refs. [5|, 16| , faceting of the surface due to strong anisotropy opposes the tendency 
of the film to dewet. Another reason is "geometrical", meaning that a planar surface of the 
as-deposited solid film may feature local defects of arbitrary shape protruding arbitrarily 
deep into the film (i.e., the pinholes). Since the attractive substrate potential decreases with 
the film height, its influence is stronger on deep pinholes, which therefore dewet faster. In 
contrast to shallow pinholes the morphology of the tip is often different from the morphology 
of other parts of the surface, i.e. the surface away from the tip may undergo formation of a 
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hill-and- valley structure due to faceting [5|, |6[ . 

These and other differences as well as importance to technologies such as the design 
and manufacture of solid thin-film devices, make dewetting of solid films a process worth 
studying. In Refs. j^, [g] analytical and computational studies are performed of the two- 
dimensional PDE-based model, which incorporates the two-layer wettin g p otential with the 
exponential decay. Previously, Golovin et al. and other authors studied similar 

models in the context of quantum dots self-assembly. Note that the two-layer potential 
model is appropriate for ultrathin solid films, while for thicker films the van der Waals 

13 1 . In this paper the model of Refs. IS 



potential has been shown to be important 



extended to the case of the two-layer wetting potential with a variable-rate algebraic decay, 
and comparisons of the two situations are performed. The models are studied using the 
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linear stability analysis, as well as the computations of the arbitrary deviation/slope surface 
dynamics. 



II. PROBLEM FORMULATION 

The governing equation for the free one-dimensional (ID) surface z = h{x,t), evolving 
by surface diffusion, has the form 

"^''^ ^;i + /4)-"^t^V (1) 



kT dx V dx^ 

where h is the height of the film above the substrate, Q is the atomic volume, D the 
adatoms diffusivity, u the adatoms surface density, k the Boltzmann constant, T the absolute 
temperature, and , = ,<«> + the surface ehemieal potential. Here ,s the regular 
contribution due to the surface mean curvature k [ij]. Also (1 + h"^) = cos^, where 9 is 
the angle that the unit surface normal makes with the [01] crystalline direction, along which 
is the 2;-axis. (The x-axis is along the [10] direction.) Thus 6 measures the orientation of 
the surface with respect to the underlying crystal structure. Note throughout the paper the 
subscripts x, t, s, u and 9 denote differentiation. 
The wetting chemical potential 



m'"' = sj(i + /4)-"'||. (2) 



where 7 is the height-dependent surface energy of the film-substrate interface. In the two- 
layer exponential wetting model [ll| 

7(/i,^) =7(/)(0) + (75-7(/)(^))exp(-V^), h>0. (3) 



In the two-layer algebraic wetting model 15| 



7iKe) = l{^^f\e)+^s)+U7^^\e)-^s)fih/i), lim/(V^) = l, lim f{h/i) = -l. 

(4) 

Here 75 = const, is the surface energy density of the substrate in the absence of the film, 
and (. is the characteristic wetting length. 7*^-^^ (6') is the energy of the film surface, assumed 
strongly anisotropic. In the exponential model 7(/;., 6) 'y^^\6) ash ^ 00, and 'y{h, 9) — > 75 
as /i ^ 0. In the algebraic model f{h/i) is such that (i) the correct surface energies, 'j^^^O) 
and 75, are recovered as h ^ ±00, and (ii) approach to the limiting value -|-1 as /i — ^> 00 is 



an algebraic power. (Of course, negative film height has no physical meaning, thus formally 
in the substrate domain h must be replaced by z in Eq. (j4]).) The suitable generic form is 



fih/i) 
which has the expansion 
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f(h/i) = i--{h/ir"^+ 



as h 
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Note that in the limit h ^ the exponential and the algebraic models give 7 = 75 and 
7 = + 7s) /2, respectively. These results follow from the 'one-sided' ('two-sided') 

nature of the the corresponding boundary layer models for the smooth transition in surface 
energy above (across) the substrate surface z = 0, over a small length scale ^. 
'j^^^O) is taken in the form 



7(^)(^) = 7o(l + cos 46) + "-k' = ^^{9) + "-k' 
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where 70 is the mean value of the film surface energy in the absence of the substrate potential 
(equivalently, the surface energy of a very thick film), determines the degree of anisotropy, 
and 6 is the small non-negative regularization parameter having units of energy. The 5-term 
in Eq. ([71) makes the evolution equation ([1]) mathematically well-posed for strong anisotropy 
17| - 23|. (The anisotropy is weak when < < 1/15 and strong when e-^ > 1/15. 5 = in 



the former case. In the latter case the polar plot of 7'^-^) (6') has cusps at the orientations that 



are missing from the equilibrium Wulff shape and the surface stiffness 7''-^^ + 7^^'' is negative 



at these orientations [2J, |25|. Thus the evolution equation is ill-posed unless regularized 
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[17|,ll8|.) The form assumes that the surface energy is maximum in the [01] direction. 
With the regularization in place, the curvature contribution to the chemical potential has 
the standard form 



(k) 



(7 + lee)K - ( y + 



(8) 



where k, = —hxxi^ + h^) s is the arclength along the surface [d/ds = {cos6)d/dx] and 
the expressions for ^{h, 0) read 



Exponential model : 7(/i, 6) = ■jp{9) + {■js — Ipi^)) exp {—h/£), 



(9) 



Algebraic model : -f{h, 6) = ^ {^^{9) + 75) + ^ (7p(^) - 7s) f{h/e), (10) 



with 7p(6') stated in Eq. ©. By using Eqs. and ([TO]) instead of Eqs. and (g]) 
we disregard the contribution of the wetting terms (exponential or inverse tangent) to the 
regularization in Eq. (jS]). Similarly, by using Eqs. and (ITU]) in Eq. ([2]), we disregard the 
contribution of the regularization term 6k'^/2 to (See Refs. s], Q] for the justification 

of this approach.) 

Using the height of the planar unperturbed film, ho, as the length scale, the nondimen- 
sional expressions for the chemical potentials read: 
Exponential model: 

/iW = (^%{e) + ^yi-exp{-h/r))K + Texp{-h/r)K-A(^^ + Kssy (Ha) 

= (g)_r) "^P(-V0 (^^^^ 

r 

Algebraic model: 

/i^'^^ = I (ipiO) + ^) (1 + / (h/r)) + ^ (1 - / (h/r)) " A (y + /..s) , (12a) 
/i(-) = ^(%W-r)f cos^, (12b) 

where 
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f{h/r) = — arctan 



1,3,5,..., 7p(^) = l + eT,cos4^. (13) 



r ^ 

Also, r = £//io is the ratio of the characteristic wetting length to the unperturbed film 
height, r = 75/70 is the ratio of the mean surface energy of the film to the substrate surface 
energy, and A = 5/ (79/19) is the non-dimensional regularization parameter. Figures 1 and 2 
show fi^^^ for both models. Note that for the algebraic model yu'-"'^ ~ for /i ^ 1, as 

follows from Eqs. (J2]), (]1]) and ([S]). 

Now using hl/D as the time scale, and the small-slope expansion in powers of e = 
\d/dx\ <^ 1, the asymptotic nondimensional evolution equation ([1]) reads: 

ht = (Pi^) - APf + P^) , (14) 



where B = VL^v^q/ [kTK^) is the Mullins coefficient and Pj?^ = —hr^^xxx- 

(1) 



In the exponential model, the terms P^ and Pw read: 

Pj:^^ = AiK^^ + A2hl.^K + A^Kxxhl - exp {-h/r) [{V + Ai) h^^^ + o.t] , (15) 
P^ = ^"^P [as/ixx/ix (1 - 5hl) + r-^h, (ai + (2a3 - 3ai) hi + a^ht)] , (16) 
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FIG. 1: The reduced nondimensional wetting chemical potential n^'^'^ = exp (— /i/r)/2r. This 
formula is obtained when 7p is taken isotropic, hx is taken zero and T = 0.5 in Eq. (jllb|) . Solid 
line: r = 0.1; dash line: r = 1. 




h 

FIG. 2: The reduced nondimensional wetting chemical potential fi^^^ = {l/4:)df /dh, where / is 
given by Eq. (jl3p . This formula is obtained when 7p is taken isotropic, hx is taken zero and T = 0.5 
in Eq. (jl2b|) . Also r = 0.1. Solid line: m = 1; dash line: m = 3; dash-dot line: m = 5. 

where Ai = ISe..^ — 1, A2 = 3 — 285e^, A3 = 2 — 150e^, ai = F — 1 — e^, 02 = F — 1 — 17e^, = 
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r — 1 + 3e^, 04 = r — 1 — 25e^. Notation o.t. (meaning other terms) in the second part of 
Eq. f|T5|) (which is proportional to the exponent and which stems from wetting interaction), 
and in the following Eqs. (1171) . (ITS]) stands for many omitted terms that do not contribute 
to linear stability. Note that the non-negative Ai signals that the surface energy anisotropy 
is strong. 

In the algebraic model the terms Pi^^ and read: 

^i^^ = 2vr (1 Hh/rr-) + + "^^^/^^''"^ [^'/^^"^ ~ 

vrae (l + (/i/r)^™)] + o.t., (17) 

= I^ai , [1 - m + (1 + m){h/r)'^] + o.t, (18) 

^ /i2 (1 + (Vr)'") 

where 05 = F — 1 + 15e^ and ag = F + 1 — 15e^. 

For computations of the surface evolution (Section IV) we use the parametric equations 

and the marker particle method. The ID surface is specified as T{x{u,t), z{u,t)), where u 

is the parameter, x and z represent the coordinates of a marker particle on a surface, which 



are governed by two coupled nondimensional PDEs 26|- 29|] : 



xt = V-Zu, (19a) 
9 

zt = -V-Xu. (19b) 
9 

Here V = B (^fii^^ + fiss ^^ is the normal velocity of the surface, and g = ds/du = \/x\\ z\ 
is the metric function. It can be easily shown that Eqs. (fT9|l are equivalent to dimensionless 
Eq. ([1]) when the surface is non-overhanging (a graph oih = hix) at all times). (Note that 
in this case u = x and d/ds = g^^d/du = (1 + hD^^^'^d / dx.) When the surface develops 
steep slope, the accurate computation using Eq. ([1]) requires a fine grid, and when the 
surface overhangs, Eq. ([T]) does not make sense. Eqs. (1191) and the marker particle method 
allow to circumvent these problems, and thus this combination is preferred for computation 
of evolving general surfaces. 

III. LINEAR STABILITY ANALYSIS OF THE PLANAR SURFACE 

A. Exponential Model 
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We assume strong anisotropy and linearize Eq. (fT^ about the equilibrium h = 1. For 
the perturbation C,{x,t) we obtain 

^t = B {Ai^^^^^ + I^^xxxxxx + exp (-1/r) [r"^ai^^^ - (r + Ai)^xxxx\ ) • (20) 

Taking ^ = e*'^^"'"'^* gives 

uj{k) = B [(Ai - exp (-1/r) (F + Ai)) - AA;^ - exp {-I / r)r-'^ aik'^] . (21) 

Note that taking the limit as r ^ in Eq. fl^T]) recovers the dispersion relation in the 
absence of wetting interaction with the substrate, uj{k) = B [Ai/c^ — AA;^]. It follows from 
Eq. (pT!) that the equilibrium surface is unstable {^^{k) > 0) to perturbations with the 
wavenumbers < k < kc, where 

kl = (2A)-MAi -exp(-l/r)(r + Ai) + 

((Ai-exp(-l/r)(r + Ai))^-4Aexp(-l/r)r-V)^^1 • (22) 

Note that the radical at the right-hand side of Eq. fl2^ always exists when ai < 0, which 
turns out to be the necessary condition for a nonwetting film |5|]. Fig. 3 shows the sketch of 
uj{k). It is interesting that the /c^-term in Eq. (pT!) coming from (Eq. (fT6l) ) makes the 
film less stable for ai < 0, but the wetting potential contribution to P^^^ in Eq. f|T5|) makes 
the surface more stable since the corresponding /c^-term is negative in Eq. ([21]) • Clearly, 
due to negative exponent this stabilizing influence is small when r is small. 

B. Algebraic Model 

Eq. gives 



u{k) = B 



2 {F + Ai + asr-^™} Aq - Trag (1 + r-^") _ 



27r (1 + r-2m 
m [1 — m + (1 + m)r~^"^^ 



-aik 



(23) 



where Aq = arctan [(1/r)™]. 

The cut-off wavenumber is compared in Fig. 4 for both models and the three values of 
m. For m = 1, kc tends to zero asymptotically, while for m = 3, 5 it becomes zero at r 
slightly larger than one. Thus the m = 1 case is qualitatively similar to the exponential 
model. Comparing the 171 = 1 case to the exponential model, it can be seen that for r < 0.5 
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FIG. 3: Sketch of the hnear growth rate uj{k). Perturbations with wavenumbers < k < kc are 
unstable and may grow nonhnearly until the film ruptures. 

the surface is more stable in the former case, and less stable for r > 0.5. Comparing the 
m = 3, 5 cases to the exponential model, it can be seen that for small values or r the 
interval of instability is the same for both models, for intermediate values of r the interval 
is larger for the algebraic model, and for r >~ 1 the surface governed by the algebraic 
model is absolutely stable, while there is still a narrow interval of long-wave instability in 
the exponential model. 



IV. NUMERICAL RESULTS FOR THE LARGE- AMPLITUDE INITIAL 
DEFORMATION (PINHOLE DEFECT) 

In this section the parameters are chosen as follows: r = 0.1, F = 0.5, = 1/12, A = 
0.005. Following the method of lines approach, Eqs. fll9p are discretized in the parameter 
u using second-order finite differences and the time-stepping is performed by the implicit 



Runge-Kutta solver RADAU 



301 ]. Initially u = x, but periodically (usually after every few 



tens of the time steps) the surface is reparametrized so that u becomes the arclength, and 
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FIG. 4: Plots of the square of the critical wavenumber vs. r. T = 0.5, = 1/12, A = 0.005. Solid 
line: exponential model; dash line: algebraic model with m = 1; dash-dot line: algebraic model 
with m = 3; dot line: algebraic model with m = 5. (Abrupt termination of the dash-dot and dot 
lines is the artifact of the plotting software. We confirmed that these lines continue to intersection 
with the r-axis.) 

the positions of the marker particles are recomputed accordingly. This prevents marker 
particles from coming too close or too far apart in the course of the surface evolution. 

We compute the dynamic morphology and its rate of evolution towards either film rupture 
or the quasiequilibrium state, which is characterized by the coarsening in time hill-and- valley 
structure at the both sides of the residual defect, which dissipates with the much slower rate 
Q]. In Ref. [sj it is shown for the exponential model that for r fixed, the outcome of the 
evolution (a rupture or a hill-and-valley structure) depends on and the initial condition, 
i.e. the width and the depth of the pinhole. As will be seen, in the algebraic model the 
outcome depends also on m, which sets the rate of change of the wetting potential. The focus 
is on the rate of the extension of the pinhole tip in the algebraic model, since the detailed 



computations for the exponential model are performed in Ref. jo], and morphologies are 
similar in both models. Also, since the parameter domain of film rupture is more narrow 
for the algebraic model, we investigate deep pinholes only. 
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The initial condition is taken as in Ref. 6|, i.e. the Gaussian curve: 



z{x, 0) = 1 — dexp 



X — 5 



w 



< a: < 10, 0<d<l. (24) 



Note that the length of the computational domain equals to ten times the unperturbed 
film height, and the defect is positioned at the center of the domain. Periodic boundary 
conditions are used. 

Figures 5 and 6 show the log-normal plots of the pinhole depth vs time, for d = 0.9 and 
w = 2, 0.15, respectively. Zm is the height of the surface at the tip of the pinhole. The 
wide pinhole dewets for m = 1 only. (See Fig. 5. Note that the exponential model predicts 
faster dewetting.) Wetting potentials with m = 3 and m = 5 result in the quasiequilibrium 
at < < 1- Quasiequilibrium means that Zm (or, equivalently, the depth) changes very 
slow or not at all, while the rest of the shape changes relatively fast. In the inset, for m = 3 
one can see the onset of the formation of the hill-and-valley structure near the endpoints 
of the domain; as has been noted in the Introduction, this does not affect the pinhole 
depth. Interestingly, here the pinhole tip is blunt at quasiequilibrium, while it is sharp in 



all examples computed for the exponential model |6|. In contrast to the wide pinhole, the 
narrow pinhole does not dewet even for m = 1 (Fig. 6) and in all three cases evolves to 
quasiequilibrium. 

It must be noted here that stable equilibrium (steady state) solutions have been numer- 
ically found in the studies of a nonlinear stress-driven morp hological instability of a solid 



film without wetting interaction, by Spencer & Meiron 3l| and by Xiang & E j32|. The 
problem under study in this paper differs from the problem studied by these authors in that 
the instability is driven not by stress but by wetting potential, and the surface energy is 
anisotropic. These instability mechanisms have different physical origins and the process 
of morphological evolution in both cases is similar but not the same. In particular, due to 
the presence of strong surface energy anisotropy the equilibrium solution, when it occurs, is 
replaced by quasiequilibrium. The latter can be viewed as the locally broken equilibrium. 
This violation of equilibrium occurs in the surface regions away from the pinhole tip. There, 
an evolving hill-and-valley structure is energetically favorable because the attraction to the 
substrate is weak. 

Finally, we note that the slight decrease of the initial depth results in the termination of 
dewetting even for wide pinholes. For instance. Figure 7 shows the case d = 0.7, w = 2. 
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FIG. 5: Kinetics (rate) data for the deep, wide pinhole {d = 0.9, w = 2). Line slope equals the rate 
of the tip evolution. Solid lines: algebraic model. Dash line: exponential model. Inset: Surface 
shapes at t = 1.8, for m = 1 (solid line) and m = 3 (dash-dot line). 



As can be seen, there is no dewetting for neither value of m. The pinhole tip is attracted 
to the substrate for a while, but then reverses the direction and will finally stabilize at a 
quasiequilibrium position. Quasiequilibrium is achieved for the m = 5 case. For comparison, 
the exponential model predicts dewetting even for the more shallow pinhole with d = 0.5 



(see Fig. 2(a) in Ref. ^\). 

To summarize, we contrasted two PDE-based models of dewetting for nonwetting ultra- 
thin single-crystal films. It remains to be seen how these models compare to experiment. 
Detailed experiments focusirig on the dynamics of a single pinhole are yet to be performed. 
(The published experiments [l|, Q] describe very briefly the initial stages of film dewetting and 
proceed to detailed study of the post-dewetting regimes, i.e. the hole widening, secondary 
instabilities and material agglomeration.) 
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FIG. 6: Kinetics data for the deep, narrow pinhole {d = 0.9, w = 0.15). Sohd hnes: algebraic 
model. Dash line: exponential model. Inset: quasiequilibrium surface shape at t = 4, for m = 1. 

I thank Brian J. Spencer for pointing out the algebraic decay model to me. 
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